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Abstract 

We study spherically-symmetric solutions in Massive Gravity generated by matter 
sources with polytropic equation of state. We concentrate in the non-perturbative regime 
where the mass term non-linearities are important, and present the main features of the 
tfs' solutions. 



1 Introduction 

Since the discovery of the late time cosmic acceleration in the 1990's [1, 2], a considerable 
amount of effort has been invested in order to understand its origin. The simplest explanation 
is the old idea of the cosmological constant which requires extreme fine-tuning in order to fit its 
value with the present matter density on one hand, and with possible microscopic mechanisms 
for its generation on the other hand [3]. This difficulty which is known as the cosmological 
constant problem motivates the search of alternative explanations for the cosmic acceleration, 
which usually are taken as one of the following two approaches: 

The first assumes some unknown form of matter with an equation of state similar to that 
of a cosmological constant: P ~ —p. This dark energy is assumed to dominate the universe at 
recent times (starting around z ~ 0.5) thus changing its expansion rate from the decelerating 
mode of matter-dominated universe to the accelerating one. Quintessence models [4, 5] are one 
popular example of this approach. 

The second approach suggests a modification of the theory of General Relativity (GR) in a 
way that produces deviations in a cosmological scale. The hope is that the modifications will 
have such an effect which will enable to explain the cosmic acceleration without invoking dark 
energy. 

These approaches may also be relevant in tackling the Dark Matter problem and the puzzle 
of coincidence, that is, that these two components, dark energy and dark matter have a similar 
weight in the present cosmological epoch. 

Massive Gravity or Massive General Relativity (MGR) [6] is a modification of GR originated 
from the very natural question: what is the way to give the graviton a mass (or a finite range)? 
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This question was partially answered back in the 1930's by Fierz and Pauli [7, 8] who wrote 
down a theory of a massive spin 2 field in Minkowski background. However, the generalization 
to the full self-gravitating case was proved to be notoriously difficult. Some of the hurdles were: 

• the van Dam-Veltman-Zakharov (vDVZ) discontinuity [9, 10] according to which gravi- 
tational attraction in MGR does not reproduce that of GR at the massless graviton limit. 

• the discovery by Boulware and Deser [11, 12] (BD) that a generic extension of this theory 
to curved backgrounds contains a ghost degree of freedom in addition to the five of the 
massive spin 2 graviton. 

A possible way out of the vDVZ discontinuity is known as the Vainshtein mechanism [13]. 
Vainshtein realized that in such a theory one cannot rely on weak field results (as vDVZ did), 
since the non-linearities become stronger as the mass (say m) of the graviton decreases. This 
phenomenon is quantified by the Vainshtein radius associated with a localized massive source 
of mass M. This new length-scale ry determines the range beyond which linear approximations 
are justified. The fact that ry increases indefinitely as m decreases (see e.g. Eq. (3.6) below) 
implies immediately that as m — )■ 0, the linear theory cannot be trusted anywhere. This opens 
the possibility that nonlinear effects will cure the vDVZ discontinuity as is indeed verified by 
an explicit analysis [14]. 

However, due to the BD ghost, MGR was considered for a long time as inconsistent and/or 
unphysical. Only very recently a solution was found [15, 16] and the theory became acceptable 
(see also [6] and references therein). 

It is quite obvious that a graviton mass term in a generally covariant theory of gravity 
cannot contain derivatives of g"^;, on one hand, and cannot make use of the only two possible 
"derivative-free" quantities tr{gii) and det{gnv) on the other hand. The conclusion is that a 
mass term cannot be based on the metric tensor alone. 

A direct way to overcome this is to introduce an additional rank-2 tensor, say Hf^^i, which 
can be used to construct together with g^y non-trivial scalars [15, 16, 17, 18]. This H^y may 
be regarded as an auxiliary metric and be parametrized by four Stxickelberg fields such 
that H^jj = Tjabdfj^Lp'^dyip^ where rjah is the Minkowski metric. Next define its square root 
by J'^Jx = H'^ and the tensor K'j^ = " J^- It was shown [15, 16] that there are only 
three combinations, quadratic, cubic and quartic that can eliminate the BD ghost in curved 
background and eventually yield a ghost-free theory. These contributions are written in terms 
of various traces of this tensor, namely /C„ = tr{K^), or actually the following combinations: 



where (T2, o"3 and (T4 are arbitrary dimensionless coefficients and the last expression uses the 
eigenvalues ka of K^. Actually, one of these (T^s can be absorbed by a redefinition of m, but 
we keep all three in order to allow each of them to vanish. 

The MGR action is thus given by the following GR contribution (k = BttG) with the 
additional three-parameter mass term (and matter Lagrangian): 



= (/Ci)2 - /C2 , >t3 = (/Ci)=^ - 3/C1/C2 + 2/C3 , 
Ma = {ICif - 6(/Ci)2/C2 + 3(/C2)2 + 8/C1/C3 - 6/C4 



(1.1) 



such that the most general gravitational mass term M is 




(^^{kokik2 + kokiks + /so^2^3 + kik2k^) + 0-4/^0^1^2^3] (1-2) 




(1.3) 
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This form of the mass term has the effect of ehminating the sixth propagating mode that 
typicahy exists in such a theory, leaving us with the five required for a massive spin 2 field. 

Although the original motivation for studying MGR was pure physics curiosity, it fits well 
to the arena of the dark energy and dark matter problems. A massive graviton is bound to 
introduce large distance effects into gravity and may have its say on the DM problem as well. 

MGR is presently in an initial state with respect to these issues, due to the harsh theoretical 
obstacles that were mentioned above. Only recently some results about Cosmology within the 
MGR framework have been published [19, 20, 21, 22]. 

As other modified theories of gravity, MGR not only introduces deviations from GR at large 
distances, but also modifies the gravitational fields around compact sources. These modifica- 
tions must be consistent with the stringent observational solar system tests, but still may be 
used in order to differentiate between the various modified gravity theories. 

A natural way to go is analyzing spherically-symmetric solution of MGR. Some preliminary 
work in this direction has been done recently [14, 24, 25, 26, 27] and here we take this direction 
and study localized self-gravitating solutions in MGR where the matter source is a perfect fluid 
with a polytropic equation of state. 

Observations already limit the graviton mass to a tiny value of 10"^^ eV at most [20, 23]. 
However, in order to understand the impact of the graviton mass term we will concentrate in 
the non-perturbative regime where the mass term non-linearities are important. This means 
that we will deal with structures whose sizes are of the order of the graviton range. 

2 Spherically-Symmetric Solutions 

For spherical symmetry we use the most general diagonal line element 

ds^ = a^{r)df - v'^{r)d? - r^c^{r){de^ + sm'^{e)dip'^) (2.1) 

and the unitary gauge where 

Vabdt,(l)''d^(t)''dxi^dx'' = dt^ - d^ - P{d9^ + sin^{e)d^'^). (2.2) 

It is well-known [14] that this way leads to asymptotically flat solutions which is the branch 
we are interested in. Note that we cannot redefine the radial coordinate to get rid of c(r) (as 
is done in GR) without breaking the Lorenz symmetry of the unitary gauge. We will have 
therefore to solve here for three components of the metric tensor. 

In these coordinates the tensor will be diagonal with components 

A;o = l-l/a, ^1 = 1-1/^;, ^2 = fcg = 1 - 1/c. (2.3) 

We can however perform a coordinate transformation r(r) such that the new line element will 
be 

ds"^ = a{rfde - dr'^/h{rf - r'^{de'^ + sin'^ {e)dip'^) (2.4) 
so the components become now 

fco = l-l/a, ki = l-b{r/cy, /c2 = /cs = 1 - 1/c. (2.5) 

The field equations obtained from the action (1.3) have the general form 

G^^ + + kT^^ = (2.6) 
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where S^^ is the contribution from the mass term. In the coordinate system where the metric 
is given by Eq. (2.4) the components of Einstein tensor are: 

_ 2bb' b^-l 

— "I — — 

Gl = ^ + ^_^ (2.7) 
ra 



^ ^ a ao ra ro 



and the mass term contribution is 



So = -/^3 + — --2+ M2- — + -2 -2 



c \ c / c- 

Sr = -fJ-a H T + ^U's (2-8) 

c a \ c ) 

^ \ac/ac\ \a c J ac J & 

where we denote /ii = vn?{o2 + 2a"3 + (74), /i2 = ?n^(3a"2 + Saa + (T4), //a = m?{6a2 + 4(73 + ^4)) 
/xo = m^(o-3 + 0-4). 

The components of the matter energy-momentum tensor may be added to the field equations 
as = diag{p{r), —Pr{r), —P±{r), —P±{r)). These components satisfy a continuity equation 

P; + {p + Pr)- + = (2.9) 

a r 

which is equivalent to the matter field equation if the matter is described by a Lagrangian. 

Otherwise, an additional equation of state is required. 

Since the matter satisfies the continuity equation separately, the mass term S^y should be 
also separately conserved. Explicitly we get 

2mi , ra' , 2(6-1) ( (^^^\^^'^\ a (0 m\ 

which is of course not an independent equation but can be useful for the analysis. 

Note the four (or actually one) functions which appear repeatedly in the equations: 

/i(a,c)=//i-Mi-if- + -V^ '^ = 2'3 (2.11) 

\a c I ac 



or more explicitly: 



/2(a,c) =//2-Mi { - + -) +— ; 52(c) = /2(c,c) 



ac 



/3(a,c)=Ai3-M2 (- + -)+— ; 53(c) = /3(c,c) (2.12) 

c ac 
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so the field equations will get the following compact form: 

2bb' ,,6(c-rc') 

262a' 62-1 , , 1 , , „ ^ 

1 n 53(c) + -52(c) - KPr=\) 

ra a 

^2 fa" a'h' a' h'\ , ,b{c-rd) 

b [ — + —r + — + -r]- Ma, c + /2 a, ^ ^ - kP^ = 

\ a a6 ra r6/ cJ 

52(c)— + A_^/2(a,c) = (2.13) 

and in expanded form: 

266' 6^-1 2ii2 Ml ( iJiQ\h{c-rc') 



H 5 /X3 H T + hW2 h -T 5 h K/9 = 

c \ c &• J & 

262a' 62-1 2/^2 m ^ 1 «o\ „ . 

^ n M3 -I 9 + - M2 1" -9 ~ = 

ra c a \ c ) 

,,/a" a'6' a' 6'\ /I l\ Mi 

6" — + -^ + — + -r -M3 + M2 - + - - — + 

a ao ra ro J \a c J ac 

(\ \\ uo\ b(c-rc') 

/X2 - Ml - + - + — - kP± = 

\ a c / ac / 



2mi Mo V«'^ 2(6-1) / l^l^lA^-^o^ . 



We may eliminate 6(r) using the last of Eq. (2.14) or alternatively, using also the second 
equation of (2.14), express a(r) in terms of 6(r), c(r) and Prir) and solve for these three 
variables. This leads to an algebraic expression 

a(r) =y(r,6,c,P,)/Z(r,6,c,P,) (2.15) 

where Y and Z are polynomials in their arguments. The equations (2.14) then reduce to a 
system of two first order equations for 6(r) and c(r). This is the approach that we followed 
to construct solutions. There exist, to our knowledge, no explicit solutions to the system 
of equations above (even for the vacuum equations), so the way to proceed was to turn to 
numerical methods and to make simple assumptions about the matter source. 

For the vacuum solutions we have however the following asymptotic behavior (for mr >> 1) 
[6, 27]: 

, 4GMe-'"'' , , 2GM e-'^Ul + mr) 

a{r) = 1 , 6(r) = 1 — 



r 3 

2GMe-"^'' / 1 1 

l + _ + 

3r \ mr m^r 



<r) = 1 + ( 1 + — + 3 ) (2-16) 



where the integration constants are expressed in terms of the ( "gravitational" ) mass M of the 
source obtained from a weak field correspondence. This will be our way to identify the mass 
of the solutions we will find. 



3 Polytrope Solutions 

We will concentrate in this work in localized asymptotically flat solutions produced by a perfect 
fluid. Asymptotic flatness is possible since a necessary condition for (a(r) 1, 6(r) 1, 
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c(r) 1) is 

/Lto - 3Aii + 3//2 - yW3 = (3.1) 

which is satisfied identicaUy. 

The simplest assumption of the nature of the perfect fluid is to use a polytropic equation 
of state that we write as 

P^ = P, = P = P,{p/P,y (3.2) 

where and 7 > 1 are parameters that characterize the matter. In order to study solutions, 
we first change to dimensionless variables. The direct way to do that is to use i = l/V kP* as 
a length scale to define x = r/£ and rh = m£ along with the obvious p = p/P*, P = P/P* and 
also fi{a,c) = fi{a,c)/m? and gi{c) = gi{c)/m?, which are still given by (2.11)-(2.12) with the 
replacements p^i ^ (li = pi/m?. The field equations become a three-parameter system whose 
solutions (a(a;), p(x) etc..) are determined by the central value of the density p(0). Typically, 
any value of ^(0) (in a certain domain) corresponds to a localized solution whose (coordinate) 
radius is determined by p{xs) = 0, that is Vg = ixg. 

However, we may turn this system into a boundary value problem which is advantageous 
from the point of view of numerical solution of this coupled system of ordinary differential 
equations. This can be done if we use as a dimensionless radial variable z = r/rg and define 
also a dimensionless radial extension e = 2{'mrs)'^- Pretending that we know in advance, the 
surface of the localized solution is now at 2; = 2;^ = 1 and we may write the field equations in 
the following way (using ' = d/dz): 

2bb' h'^ -I e /_ , , _ , Xc - zd) 
— + ^ - 2 ( ^3(c) - -92{c)^^ 



26V e( 1 . _^ 
+ — 2 o 53(c) - -52(c) - ropT 



za ' z"^ 2 V^"'"'^ ' a 

,ofa" a'b' a' b' \ e/-, , ,b(c-zc')\ ^ ^ 

b^- + -r + — + -r)-7^(h{a,c)-f2{a,c) ^ ^ ' )-wfP = (3.3) 
\ a ab za zb / 2 \ & ) 

, ^zo! 2(6- 1) - , , ^ 
52(c)— + ^-^/2(a,c) = 

7P +(p+(p)'-^)- = 0. 

a 

An additional parameter, w = kP^t]. = {rs/tf appears in the equations, but now wc expect 
that for each (ro.e./Jj) there will typically be a single solution {a{z), p{z) etc..) that satisfies 
the simple boundary conditions p'{0) = and p(l) = 0. 

The solutions can be characterized in terms of a few physical parameters, namely the mass 
M and size of the solutions. We may define also the "inertial mass" Mj and the physical 
radius given by 

r r - £ - f^" dr 

GMi = G d^x^/\i\T^ = -^Miw = -Miw^/^ , r^^ys = ttt = iz^hys^^/^ (3.4) 

Jr<rs II Jq 0[r} 

where we use the dimensionless quantities 

r z^a{z)p{z) dz ,„ 

Occasionally we will use the Schwarzschild radius th = 2GM in order to present the mass and 
to compare with GR results. We will calculate also the Vainshtein radius given for this theory 

by [6] 

rv = (87rGM/m2)V3 (3.6) 
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where we use the "gravitational mass" M which is anyhow quite close to Mj. The way to 
obtain M is by matching our solutions to the asymptotic vacuum behavior given in (2.16). 

We will see that for the parameter space that we study, the Vainshtein radius is not large 
enough in order for the Vainshtein mechanism to be realized. 

4 Numerical Techniques and Vacuum Solutions 

As mentioned already, it seems that no explicit solution to the above system of equations can 
be obtained, therefore we have used a numerical technique to construct the solutions. The 
radius Zg and the parameters appearing in the mass term have, of course, to be fixed a priori. 
Without loosing generality, we can fix the scale and set Zs = 1. We need also to specify the mass 
parameters and we chose the representative values used already by Gruzinov and Mirbabayi 
[25] /ii = 2 /x2 = 3 and /X3 = 5. Note however that the solutions of Ref. [25] correspond to a 
source of incompressible fluid, with a much smaller value of e than we use here. 

The field equations were solved by employing a collocation method for boundary-value 
ordinary differential equations, equipped with an adaptive mesh selection procedure [28]. The 
most efficient way we found to construct solutions consists of solving the equations in two steps 
that we now describe : 

• First we solve the vacuum equations on an interval z £ [zg, 00] by imposing the boundary 
conditions 

c{zs) = Cs , 6(00) = 1 (4.7) 

where Cg is a constant and where the second condition ensures the solution to be asymp- 
totically flat. The values bg = b{zs) and = a{zs) of the metric fields corresponding 
to Cs can then be extracted from the numerical solution. Actually, some c<j values may 
correspond to more than a single vacuum solution, so in Fig. 1 we plot and Cs as a 
function of bg. The plot is for e = l/2, e = 2 and e = 8. Not shown in the plot are the 
GR curves, = bg and = 1. More generally, we expect such families of asymptotically 
flat solutions (labelled e.g. by bg) to exist for generic values of the mass parameters /ij. 

• For the second step, i.e. solving for < z < z^, it turns out essential to impose four 
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U , , , 1 , , , 1 , , , 1 , , , 1 , , , L 

0.0 0.2 0.4 0.6 0.8 1.0 

*. 

Figure 1: Values of (lower curves) and of Cg (upper dashed) as function of bg for the vacuum 
solution with e = 1/2, 2, 8. The GR solutions correspond to the lines = bg and Cg = 1. 
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boundary conditions for our numerical methods to work efficiently, namely 

6(0) = 1 , p(0) = po , P(^.) = , c{zs) = cs (4.8) 

where po is a positive constant. The first condition is necessary for the metric to be 
regular at the center. Two conditions on p{r) turn out to be necessary, otherwise the 
linear equation of the density function leads to the trivial solution p{r) = 0. The last 
condition is for the continuity of the field c{z) at z = Zg. Since the system consists of 
three first order equations, it can accommodate only three boundary conditions. In order 
to be able to impose (4.8), the issue consists in extending the system to four equations 
by supplementing the equation dw/dr = 0. The parameter w is therefore not an input 
in this procedure, but is a part of the solution, thus fixing the radial size of the solution. 
To proceed, we first choose a couple of values (65,05) determined in the first step and 
then solve the system with the conditions (4.8). The value pQ has to be fine tuned in 
such a way that h{zs) = bg. The fact that the equations for b{r),c{r) are first order then 
guarantee the metric functions to be smooth at z = Zg. 

5 Solutions in Massive Gravity 

The numerical integration of the equations turns out quite involved and we did not address the 
system throughout the full range of parameter values. We discuss here the solutions for the 
polytropic power 7 = 2 and two values of the parameter e which reflect the main features of the 
solutions. We consider first a case where the radial extension parameter is large: e = 8. Then 
we report the results obtained in a case where the radial parameter is intermediate, i.e. e = 2. 
Surprisingly, our numerical technique works better for large e and the pattern of solutions is 
obtained more easily in this case. As mentioned above, we set = 1 by an appropriate scaling. 

5.1 Solutions in General Relativity 

In order to test our numerical procedure and for the sake of comparison of the new results, we 
first solved the equations in the case of GR. In this case the independent fields are a{r),b{r) 




12 3 4 

P(0) 



Figure 2: Central density dependence of parameters of the GR solutions: Tphyg, rg, rjj — 2GM 



8 



and p{r) and the corresponding (first order) equations have to be solved with the following 
boundary conditions : 

6(0) = 1 ,p(r,) = , a{r,) = b{rs) (5.9) 

By the rescaling explained above will correspond to z = r/rg = 1, i.e. Zg = 1 and the 
only free parameter is zu. The solutions can then be constructed for different values of zu. 
Our numerical results show that the solutions exist for values of tz7 in a finite interval, i.e. for 
zurn < < wm (wc find zum ~ 6 and wm ~ 40) with a single solution for each w. In the 
limit w — )■ Wm, the density function approaches p{r) = and vacuum solution (i.e. Minkowski 
space-time) is recovered. In the limit w — )• zUm, we observe that the value of metric component 
500 at the center, i.e. a(0) tends to zero and the solution becomes singular. The values of 
density at the center p{0) diverges, but the mass stays finite. Fig 2 shows the behavior of the 
main parameters of the solutions as functions of the central density p{0) = p{0)/P^. 

5.2 Massive Gravity Solutions: Case e = 8 

Next we move to the main subject of this work, the MGR solutions. Our results are summarized 
in Fig. 4 where several parameters characterizing the solutions are plotted. The main structure 
is seen from Fig. 4 -left. The right hand part shows further details. 

One new feature is that the number of values of p{0) leading to a matching of the metric 
functions at the surface, i.e. to b{zs) = bs,c{zs) = Cg depend very sensitively on the values 
of e and of bg or Cg. Actually, only the small rightmost interval of Fig. 1 is realized in the 
present case: For bg < 0.845 (which correspond to Cg > 1.089), the values b{zg) computed from 
the interior solutions are always larger than bg and no continuous solution is available. In the 
interval 0.845 < 6s < 0.899 (1.058 < Cg < 1.089) there are two solutions corresponding to the 
same values of {bg,Cg), they are distinguished by the values of p{0). That is two different interior 
solutions correspond to the same exterior solution. Profiles of two such solutions for (6^ = 0.870, 
Cg = 1.075) are illustrated in Fig. 3 (left part). The two solutions have {p{0) 1.41, zu 21.2, 




Figure 3: Profiles of two solutions with the same exterior behavior corresponding to e = 8, 
bg = 0.870, Cg = 1.075: Left: The curves for p{0) = 0.26 (red) and p{0) = 1.41 (black). Added 
for comparison is the exterior a(r) = 6(r) curve of GR with the same mass and as the 
second solution: GM/i = 2.39, rg/£ = 5.96. Right: Fit of the exterior solutions (solid lines) 
by the asymptotic solutions (2.16) (dashed lines). Beyond r/rg = 2 the difference between the 
solutions cannot be resolved in the figure. 
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P(0) ' „, 

Figure 4: Characteristics of the solutions for e = 8. Left: Mass and radii of the solutions as 
functions of the density at the origin /9(0): M, r<j, ry and r^hys . The dashed curves are the 
corresponding GR ones. Right: The same parameters plus /)(0) vs hg. 

GM/l 1.85) and (p(0) w 0.26, w k, 35.5, GM/l ^ 2.39). The right-hand part of this figure 
demonstrate the very good fit of the exterior solution by the asymptotic solution (2.16) starting 
already around r/r^ = 1.5. This figure shows the domain where the actual and perturbative 
solutions are still discernable. 

The existence of two branches occurs also in GR as can be inferred from the GR plots (Fig. 
2), but it is more pronounced in MGR. At any rate, this does not violate any basic principle, 
since the two different interior solutions are actually solutions of two systems differing by the 
values of w. From a physical point of view, these two solutions have different sizes, so it is only 
when they are written in terms of the ratio r/rg, that the exterior solutions match. 

The two branches of solutions do not persist beyond bs = 0.899 (or Cg = 1.058). Indeed 
we observe that the value a(0) associated with one of the branches decreases and becomes null 
in the limit bg = 0.899. This is of course the limit p(0) — )■ oo mentioned above - see Fig. 4 
right. Therefore only one solution exists for any 0.899 < bg < 1 (or 1 < Cg < 1.058). In the 
limit {bg — )• 1, — )• 1), the matter density converges uniformly to p{r) = and the limiting 
configuration is just the vacuum. 

As for the parameter w = (rg/i)'^ defined above, we find that only a subset of the values 
of the parameter vj allowing for polytrope solutions in GR lead also to MGR solutions. Also 
there is a small interval of tu (for instance 19.8 < w < 20.4 which corresponds to p{0) ~ 4.5) 
for which two solutions exist - see the region near the local minimum of the rg/i curve in Fig. 
4 -right. In other words, in a certain domain in the parameter space two values of the central 
density correspond to the same value of rg. The plot of rg/i vs. p(0) develops small oscillations 
which cannot be resolved in Fig. 4 -left. This feature seems peculiar to MGR. As another 
feature, we observe that, in the critical limit where a(0) tends to zero, p{0) diverges while the 
inertial mass remain finite as in the GR case. 

Last we turn our attention to the Vainshtein radius ry. As is obvious from Fig. 4, its 
maximum value is around ry ~ 8 which is never significantly larger than rg. We therefore 
do not expect the Vainshtein mechanism to operate for this family of solutions. On the other 
hand, the perturbative (weak field) solution - Eq. (2.16) approximates very well the solutions 
as is demonstrated in Fig. 3 -right. 
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Figure 5: Characteristics of the solutions for e = 2. Left: Mass and radii of the solutions as 
functions of the density at the origin p{0): M, r<j, ry and rp^yg . The dashed curves are the 
corresponding GR ones. Right: The same parameters plus p{0) vs bs- 

5.3 Massive Gravity Solutions: Case e = 2 

Solutions in this case exist for a larger interval of values of bs or Cg- However, it became quite 
difficult to construct them. The reasons for this difficulties is technical and results from the 
elimination of the field a(r) in the equations as explained at the end of Sec. 2. Indeed, it turns 
out that when the parameter Cg becomes large, both numerator and denominator (denoted X 
and Y in (2.15)) develop a zero at some common value of the radial variable, say r = rg. This 
creates numerical difficulties although the resulting field a(r) is continuous and differentiable. 
We checked that this problem also occurs when b(r) is eliminated instead of a(r) from Eqs. 
(2.14). Avoiding this problem would likely require to adopt a different coordinate system but 
this direction was not taken in the present paper. Anyhow, the pattern of the solution is similar 
to the case e = 8 and the counterpart of Fig. 4 for the case e = 2 is included as Fig. 5. Solutions 
exist now in a somewhat larger interval of 0.765 < 6s < 1 where bg = 0.818 divides between the 
two branch domain and the one branch domain: for 0.765 < bg < 0.818 two interior solutions 
exist for a single exterior one, while for 0.818 <bs<l there exists a single branch. 

Here too, the Vainshtein radius ry is not large enough for turning on the Vainshtein mech- 
anism. As is obvious from Fig. 5, its maximum value is around ry ~ 9 which is never 
significantly larger than r<j. On the other hand, the perturbative (weak field) solution - Eq. 
(2.16) still approximates very well the solutions very much as demonstrated in Fig. 3 (right 
part). 

5.4 Other Cases 

The /9(0)-dependence of the mass and radial parameters are shown for these two values of e 
together with the GR curves in Fig. 6. 

We further studied the system also for smaller values of e and obtained results qualitatively 
similar to the case e = 2. In particular the Mass parameters M, Mj approach continuously their 
corresponding GR values. In all cases we find a similar behavior of a maximal mass for a certain 
central density beyond which the mases decrease and the solutions are unstable. However, it 
seems that the graviton mass term increases considerably the mass range of possible solutions. 
The Vainshtein radius increases as e decreases and we expect this tendency to continue such that 
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Figure 6: Characteristics of the solutions for e = 2, 8 together with the GR curves (dashed). 
Left: The mass parameters Mj (lower blue lines) and M (upper red lines) as a function of the 
central density. Right: The radial parameters (lower blue) and ry (upper red). 

the Vainshtein mechanism will take action and a Newtonian domain will emerge at intermediate 
distances from the source. 

Finally, we studied the equations in the case when only the quadratic mass term is present, 
i.e. £73 = £74 = 0. Preliminary results obtained in this case confirm that solutions exist roughly 
with the same pattern as in the previous cases that were studied more extensively. However, 
the construction also reveals that the numerical difficulties related to the occurrence of zeroes 
in the numerator and denominators of o(r) in Eq. (2.15) are more severe than in the case 
where the higher order terms are present. 

6 Conclusion 

This paper summarizes the first steps of analysis of spherically-symmetric solutions in a recently 
proposed ghost-free model for massive gravity (MGR), where the non-linear effects cannot be 
treated perturbatively. 

We have constructed new families of spherically-symmetric solutions in MGR generated 
by perfect fluid matter sources. The solutions are obtained non-perturbatively by solving the 
underlying non-linear equations, they are regular over the whole space-time, concentrated inside 
a sphere of radius r^, and characterized by an energy density and a pressure related by a power 
law: a polytrope. The equations are a system of coupled non-linear equations labeled by two 
physical parameters: the graviton mass (or inverse range) m and the length scale i defined 
in terms of the energy density scale hy i = l/^/SnGP^ which enters naturally into the 
equations. 

Within the range of parameters for which we have been able to construct numerically 
reliable solutions, we solved the equations along curves in parameter space of constant mr^. 
The solutions were obtained for a rather limited domain where mr^ is of order one. We believe 
that this limitation of our solutions is just a technical problem which can be overcome by 
more sophisticated solving techniques. This is of paramount importance since a more extensive 
analysis is evidently required in order to go beyond the e ~ 1 range. 

We found that the solutions exist only for values of the ratio rg/i limited to quite a small 
interval. This is the reason why the Vainsthein radius turns out to be not much larger than 
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in the region of parameter space we mapped, and the Vainshtein mechanism is not operative. 
In other words, the hmit of vanishingly smaU m is out of reach for the parameters we used. 

For a given central density p(0) the mass and size of the solutions are larger than the 
corresponding GR ones and increase with e. The Vainsthein radius on the other hand decreses 
with e pointing to the possibility that the Vainshtein mechanism turns on for small enough 
values of e. 

As mentioned already, studying this system for small enough e is a major part of this 
program which is currently in progress. Also interesting is a systematic survey of solutions for 
the full 2-dimcnsional parameter space of mass parameters o"j. It is already known that the 
case ((72 = 1, (J3 = —1, (J4 = 1) is special in that the Vainshtein mechanism does not operate 
[29] since the system is "not nonlinear enough" . It is expected that several more special points 
in this parameter space exist which may shed light on its structure and hint towards a preferred 
mass term for MGR. 

Also the value 7 = 2 of the polytropic power was chosen for convenience, and extending 
our work for 7 > 2 as well as for 1 < 7 < 2 is a self-evident direction for further study. 
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